module DQ
    implicit none
    integer(4), private, parameter :: iwp=SELECTED_REAL_KIND(15)
    contains
        
    subroutine DQCoeff(xlist, maxorder, c)
    !xlist: input, one-dimensional real list, with length n
    !maxorder : input, interger
    !c : output, differential quadrature coefficient array, with dimensions (0:maxorder, n, n)
    implicit none

    integer(4) :: maxorder, i, j, m, n
    real(kind=iwp) :: xlist(:)
    real(kind=iwp) :: c(0:,:,:)
    
    n = ubound(xlist,1)
    c = 0.0_iwp
 
    do i = 1, n
        do j = 1, n
            if (i /= j) then
                c(1,i,j) = Mp(xlist,i)/((xlist(i) - xlist(j))*Mp(xlist,j))
            end if
        end do
    end do
    
    do i = 1, n
        c(1,i,i) = -sum(c(1,i,1:i-1)) - sum(c(1, i, i+1:n))
    end do
    
    do m = 2, maxorder
        do i = 1, n
            do j = 1, n
                if (i/= j) then
                    c(m,i,j) = m * ( c(m-1,i,i)*c(1,i,j) - c(m-1,i,j) / ( xlist(i) - xlist(j) ) )
                end if
            end do
        end do
        do i = 1, n
            c(m,i,i) = -sum(c(m,i,1:i-1)) - sum(c(m,i,i+1:n))
        end do
    end do
    
    do i=1,n; c(0,i,i) = 1.0_iwp; end do
    end subroutine DQCoeff


    function Mp(xlist, i)
        implicit none

        real(kind=iwp) :: Mp, xlist(:)
        integer(4) :: n, i, j
        n = ubound(xlist,1); Mp = 1.0_iwp
        do j = 1, i-1
           Mp = Mp * ( xlist(i) - xlist(j) )
        end do
        do j = i+1, n
            Mp = Mp * ( xlist(i) - xlist(j) )
        end do
    end function Mp

    subroutine expandedCheb(n, xlist)
        ! expanded chebyshev points
        implicit none
        real(kind=iwp), parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592_iwp
        integer(kind=4) :: n, k
        real(kind=iwp) :: xlist(n), div

        if (n == 1) then
            xlist(1) = 0.0_iwp
        else
            div = cos(pi / (2.0_iwp * n))
            do k = 1, n
                xlist(k) = cos((2.0_iwp * (n - k + 1.0_iwp) - 1.0_iwp) * pi / (2.0_iwp * n)) / div
            end do
        end if

        return
    end subroutine expandedCheb

    subroutine LocalDQCoefficientsOdd(arr, L, maxorder, ret)
        !ret(0:maxorder, n, n); c(0:maxorder, L, L); arr(n)
        implicit none
        integer(kind=4) :: i, m, n, L, maxorder
        real(kind=iwp) :: arr(:), ret(0:, :, :)
        real(kind=iwp), allocatable :: c(:, :, :)

        ret = 0.0_iwp
        n = ubound(arr, 1)
        allocate(c(0:maxorder, L, L))

        if (L > n) then
            write(*,*) "Error! too large local point number L. L should be no more than total node number."
            return
        elseif (mod(L,2) == 0) then
            write(*,*) "Error! L should be odd number."
            return
        endif

        call DQCoeff(arr(1 : L), maxorder, c)
        do m = 0, maxorder
            ret(m, 1 : (L + 1)/2, 1 : L) = c(m, 1 : (L + 1)/2, :)
        end do

        do i = 1 + (L + 1)/2,  n - (L + 1)/2
            call DQCoeff(arr(i - (L - 1)/2 : i + (L - 1)/2), maxorder, c)
            do m = 0, maxorder
                ret(m, i, i - (L - 1)/2 : i + (L - 1)/2) = c(m, (L + 1)/2, :)
            enddo
        enddo

        call DQCoeff(arr(n - L + 1 : n), maxorder, c)
        do m = 0, maxorder
            ret(m, 1 + n - (L + 1)/2 : n, n - L + 1 : n) = c(m, (L + 1)/2 :, :)
        enddo
        return
    end subroutine LocalDQCoefficientsOdd

    subroutine LocalDQCoefficientsEven(arr, L, maxorder, ret) 
        !ret(0:maxorder, n, n); c(0:maxorder, L, L); arr(n)
        implicit none
        integer(kind=4) :: i, m, n, L, maxorder
        real(kind=iwp) :: arr(:), ret(0:, :, :)
        real(kind=iwp), allocatable :: c(:, :, :)

        ret = 0.0_iwp
        n = ubound(arr, 1)
        allocate(c(0:maxorder, L, L))
        if (L > n) then 
            write(*,*) "Error! too large local point number L. L should be no more than total node number."
        elseif (mod(L,2)==1) then
            write(*,*) "Error! L should be even number."
        endif

        call DQCoeff(arr(1 : L), maxorder, c)
        do m = 0, maxorder
            ret(m, 1 : 1 + L/2, 1 : L) = c(m, 1 : 1 + L/2, :)
        end do
        do i = 2 + L/2, n - L/2 + 1
            call DQCoeff(arr(i - L/2 : i + L/2 - 1), maxorder, c)
            do m = 0, maxorder
                ret(m, i, i - L/2 : i + L/2 - 1) = c(m, 1 + L/2, :)
            enddo
        enddo

        call DQCoeff(arr(n - L + 1 : n), maxorder, c)
        do m = 0, maxorder
            ret(m, n - L/2 + 2 : n, n - L + 1 : n) = c(m, L/2 + 2 :, :)
        enddo
        return
    end subroutine LocalDQCoefficientsEven

    subroutine LocalDQCoefficients(arr, L, maxorder, ret)
        !ret(0:maxorder, n, n);arr(n)
        implicit none
        integer(kind=4) :: L, maxorder
        real(kind=iwp) :: arr(:), ret(0:, :, :)

        If (mod(L,2)==0) then 
            call LocalDQCoefficientsEven(arr, L, maxorder, ret)
            return
        else
            call LocalDQCoefficientsOdd(arr, L, maxorder, ret)
            return 
        endif
    end subroutine LocalDQCoefficients
    
end module DQ